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Abstract 

In the fermion loop formnlation the contributions to the partition function naturally 
separate into topological equivalence classes with a definite sign. This separation forms the 
basis for an efficient fermion simulation algorithm using a fluctuating open fermion string. 
It guarantees sufficient tunnelling between the topological sectors, and hence provides 
a solution to the fermion sign problem affecting systems with broken supersymmetry. 
Moreover, the algorithm shows no critical slowing down even in the massless limit and 
can hence handle the massless Goldstino mode emerging in the supersymmetry broken 
phase. In this paper - the third in a series of three - we present the details of the 
simulation algorithm and demonstrate its efficiency by means of a few examples. 


1 Introduction 

The reformulation of supersymmetric quantum mechanics on the lattice in terms of bosonic 
and fermionic bonds as derived in the first paper of our series [1] provides a perfect setup for 
Monte Carlo simulations. First of all, the reduction in complexity by going from continuous 
to discrete variables is enormous. More specifically though, expressing the Grassmann fields 
in terms of fermionic bonds avoids the expensive calculation of the fermion determinant and 
allows the use of special algorithms for which critical slowing down is essentially absent [2, 3] 
and simulations are possible even in the massless limit [4]. This is of particular importance 
for systems with broken supersymmetry, since the physics of those is driven by the massless 
Goldstino mode. In the present paper - the last in a series of three - we describe in detail 
such an algorithm and demonstrate its efficiency. Since the model can be solved exactly at 
hnite lattice spacing by means of transfer matrices, as discussed in the second paper of our 
series [5], there is in principle no need for numerical simulations. Hence, the present paper 
rather constitutes a feasibility study to test the practicability and efficiency of the proposed 
simulation algorithm for the quantum mechanical system in the bond formulation. In that 
sense it also serves as a preparation for the application of the algorithm, in particular the 
fermionic part, in more complex situations, such as in supersymmetric Yang-Mills quantum 
mechanics [6], in the M = 1 Wess-Zumino model [7-9] or in the supersymmetric nonlinear 
0{N) sigma model [10]. The advantage of the application of the algorithm in the quantum 
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mechanical model presented here is of course the fact that the correctness of the algorithm 
can be crosschecked with the exact results from the transfer matrix approach, and that the 
algorithm can hence be validated in detail. 

There is another rather pedagogical reason which motivates to consider a new simulation 
algorithm for quantum mechanics in the bond formulation. Often, simple quantum mechanical 
systems such as the harmonic and anharmonic oscillator are used to introduce the path 
integral approach. Similarly, the systems also provide a pedagogical context in which various 
Monte Carlo simulation algorithms can be illustrated and discussed, see for example [11] 
for an early example. However, it turns out that the standard Metropolis algorithms and 
even more advanced algorithms such as the overrelaxation or heat bath algorithm become 
extremely inefficient towards the continuum limit. This has to do with the usual critical 
slowing down of the simulations towards that limit, and for the anharmonic oscillator also 
with the suppressed tunnelling at small lattice spacing. The algorithms presented here do 
not suffer from these deficiencies, because they eliminate critical slowing down. In addition, 
in the bond formulation the Z 2 -symmetry cj) —)• —cj) is exactly maintained for each bond 
configuration. 

Last but not least, the numerical simulations presented here serve as a test of the practi¬ 
cability of the solution of the fermion sign problem proposed in [4] and discussed further in 
the first paper of our series [Ij. The solution is based on two ingredients. Firstly, the lat¬ 
tice regulates the vanishing Witten index and therefore also the sign problem. Secondly, the 
fermion loop formulation provides a tool to handle the fluctuating sign, because it naturally 
separates the contributions to the partition function into topological equivalence classes, each 
possessing a definite sign. Nevertheless, it is a priori not clear whether the lattice artefacts 
and the statistical fluctuations can be kept under sufficient control in a practical simulation. 
The statistical fluctuations of the sign are essentially determined by the amount of tunnelling 
between the topological sectors, i.e., between the fermionic and bosonic vacuum. In order for 
the fermion update algorithm to be a true solution to the sign problem, it must guarantee a 
sufficiently efficient tunnelling rate. The results in this paper demonstrate that this is indeed 
the case. Not surprisingly, the open fermion string algorithm discussed here has proven to be 
extremely successful in the J\f = 1 Wess-Zumino model [9] in which the fermion sign problem 
is prevailing. 

Of course, supersymmetric quantum mechanics has already been simulated on the lat¬ 
tice in various setups using standard algorithms, cf. for example [12-20]. However, the bond 
formulation together with the simulation algorithm presented here brings the numerical non- 
perturbative calculations to a new, unprecedented level of accuracy. In that sense, the results 
presented here and partly in [4] serve as a benchmark against which new formulations or 
simulation algorithms can be tested. 

The present paper is organised as follows. In Section 2 we construct in detail an algorithm 
designed for updating the bosonic and fermionic bond configurations. The discussion includes 
the explicit update steps and the derivation of the corresponding acceptance ratios. Their 
evaluation requires the calculation of site weight ratios which turn out to become numerically 
unstable for large site occupation numbers. Therefore, in Section 3 we present a computational 
strategy which allows to evaluate the ratios for arbitrarily large occupation numbers. In 
Section 4, we then present the results obtained using the proposed algorithm. The simulations 
are for the same discretisation schemes and superpotentials we used in the previous two papers 
[1, 5]. Since this section is merely meant as a validation of the algorithm, the discussion of 
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the physics behind the results is kept short and we refer to the exact results in [5] for a more 
thorough discussion. 


2 Simulation algorithm 

We start our discussion from the partition function of supersymmetric quantum mechanics 
on the lattice written as a sum over all allowed, possibly constrained bond configurations 
C = {n\{x),n^[x)} in the configuration space Z, 

Z=J2Wf{C), (1) 

ccz 

where the fermion number F = 0,1 is determined by the fermionic bond configuration {n-^ (x)} 
with {x) = F, Vx, and the weight Wf{C) of a configuration is given by 

( 2 ) 

2 ^ / tC 

Here, x denotes the sites of the lattice and i labels the various types of bosonic bonds 
with i G {j ^ k\j,k G ¥\}. The corresponding bosonic bond weights are denoted by Wi and 
n\{x) G No is the occupation number of the bond bi connecting the sites x and x +1. The site 
weight Qf depends on the site occupation number, i.e., the total number of bosonic bonds 
connected to site x, 

N{x) = X] (-^ ■ + k • n5^fc(x - 1)) (3) 

j,k 

and is given by 

/ OO 

d(t) . (4) 

-OO 

In Section 3 we will discuss in detail the computational strategy necessary to reliably evaluate 
ratios of these integrals for arbitrary and possibly large site occupation numbers. The type of 
bonds bi, the weights wi as well as the potential V{4') and the monomer term M(4>) in eq.(4) 
depend on the specifics of the chosen discretisation and the superpotential P{4>). We refer to 
the appendix of our first paper [1] for a compilation of the discretisations and superpotentials 
considered in our series. 

As mentioned above, the bond configurations C = {n^{x), nk(x)} are possibly constrained. 
In particular we have the local fermionic constraints 

{x — 1) = {x) (5) 

while the local bosonic constraints 

N{x) = 0 mod 2 (6) 

may or may not be present depending on the bosonic symmetries of the system. 

The challenge of updating constrained bond configurations lies precisely in the difficulty 
to maintain the constraints while moving efficiently through the configuration space Z. In 
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[21] Prokof’ev and Svistunov proposed to extend the constrained bosonic bond configuration 
space by introducing local sources which explicitly violate the constraints. The so-called 
worm algorithm then probes the extended configuration space by moving the local violations 
around the lattice, thereby sampling directly the bosonic correlation function corresponding 
to the sources introduced. The contact with the original configuration space Z is established 
when the violations annihilate each other, e.g. when moving to the same site on the lattice, 
such that the bond configuration fulfils again all constraints. 

In [2] the idea has been extended to fermionic systems expressed in terms of fermionic 
bonds. The fermionic constraint in eq.(5) allows only either an empty or a completely filled 
fermion bond configuration. The difficulty for the direct application of the worm idea to the 
fermionic system lies in the fact that the introduction of the fermionic source term is 

incompatible with the presence of the fermion loop at site x. A simple solution is to allow the 
unphysical situation of the site x being occupied by a propagating fermion and two additional 
sources. Such a configuration violates the Pauli exclusion principle and does not contribute 
to any physical observable. In the Grassmann path integral such a configurations indeed 
vanishes trivially. 

In order to be more explicit it is necessary to introduce the bond configuration spaces of 
bosonic and fermionic two-point correlation functions, Qp and , respectively, following the 
notation in our first paper [1]. Bond configurations in Qp contribute to the non-normalised 
bosonic two-point function according to 


while the configurations in contribute to the non-normalised fermionic two-point function 
as 


aHxi - X2) = ((</>!,' 01 ,)) = 


E n 

ccgf ix&r 


Qi{N{x)) 

Qo{N{x)) 


• Wo{C ), 


( 8 ) 


where J- denotes the set of lattice sites belonging to the open fermion string associated with the 
fermionic correlation function. The key point of the bosonic and fermionic updating algorithm 
is that the bond configurations for g'^(O), (0) and Zp have identical bond elements. As a 

consequence, statistics for and Z can be accumulated in the same Monte Carlo process. 
If the bosonic constraints in eq.(6) are not present, e.g. for superpotentials with broken 
supersymmetry, the equivalence of bond configurations even extends to g^p{x), i.e., Zp = Qp. 
The movements from one configuration space to the other are induced by introducing or 
removing bosonic or fermionic sources according to the scheme given in figure 8 of our first 
paper [1]. For convenience we reproduce it here in figure 1. 

In the following we will now discuss in detail the various updating steps which establish 
explicitly the connection between the bond configuration spaces Qf, Zp,Qp and in addition 
move the system within and Qp. The moves are generated by a Monte Carlo process with 
probabilities given by the weights of the configurations in eqs.(2), (7) and (8). In particular, 
we derive the transition probabilities Px{C —?■ C') for the transition X from bond configuration 
C to C', which is then accepted by the usual Metropolis prescription 


PacciC ^ C') = min{l, Px{C ^ C')}. 


(9) 
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Figure 1: Schematic representation of the configuration spaces. The configuration space = Qq — Q{ 
mediates between the bosonic and the fermionic sector. By the symbols © and 0, we denote the addition and 
removal of the corresponding source and sink field variables, respectively. 


In order to simplify the discussion we select the update from Zp to Q^p or with equal prob¬ 
ability which is balanced by corresponding proposal probabilities to select between moving in 
and Q^p or returning to Zp. 


2.1 Updating the fermionic bond confignration 

Here we discuss the various update steps which moves the system within and relate the 
bond configurations spaces Zq and Zi via . 

Moves within are induced by shifting ijj by one lattice spacing from site x to site x -|- 1, 
and vice versa, while keeping the other source if) hxed. Such an update step is graphically 
illustrated in figure 2 and is called ‘shift’ update step. A shift in forward direction, x —>■ x-|-l, 


X + 1 ^ X 


( 


■0 

X-*--- 

X \ 

— X ^ X + 1 

0 / 

o X- 

X + 1 


Figure 2: Fermionic bond configuration update algorithm. Graphical representation of the ‘shift’ update step 
X —>■ X + 1 in forward direction for an open fermion string configuration. It is balanced with the shift update 
step X + 1 —>■ X in backward direction. The bosonic background bond configuration is not drawn. 


automatically involves the removal of the fermionic bond b-^ (x), whereas a shift in backward 
direction, x -|- 1 —)• x, requires the addition of a new fermionic bond {x). Both directions 
are proposed with equal probability 1/2 and are hence balanced against each other as long 
as the new site does not coincide with the position of the source 0. The formula in eq.(8) 
provides us with the acceptance ratios 


^ ' <3o(ivp))’ 

(10) 

QiyN{x)) 

(11) 
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Now let us consider the case when we propose to shift the source V’ forward to the site x 
where the source V' is present, as depicted in the upper half of figure 3. 

t/i -i/) 


X--© 



o o o o 


Figure 3: Fermionic bond confignration update algorithm. Graphical representation of the ‘shift’ update step 
x—^x—loix — l—^x, respectively, and the ‘put/remove’ update step 0 —> a; and 0 ^ x. The sources are 
marked with a O for ‘ip and a x for ip. The bosonic background bond configuration is not drawn. 


The forward shift update step x — 1 —)■ x is balanced with the backward shift update 
step X — )• X — 1. This backward shift, however, is proposed with probability 1 instead of 
probability 1/2 since the shift of ip from x —x + 1 would involve the creation of an open 
fermion string around the entire lattice. The asymmetry in the proposition probabilities is 
balanced by the choice of the probability = 1/2 to remove the sources ipip, such that we 
find the acceptance ratio for a shift in backward direction x —)■ x — 1 to be the same as in 
eq.(lO), namely 


Psh{x X - 1) 


Qi(iV(x-l)) 

Qo(iV(T-l))' 


( 12 ) 


The shift step is balanced with the corresponding one in forward direction with acceptance 
ratio as given in eq.(ll). 


Psh{x 


Qo(jV(x-l)) 
^ Qi(iV(x-l)) 


(13) 


The step from to Zq and vice versa is induced by introducing or removing a pair of 
fermionic sources ipip at site x, respectively. It is called ‘put/remove’ update step and is 
graphically illustrated in the lower half of figure 3. The removal of the fermionic sources 
is suggested with probability p^m = 1/2 and is balanced on one side by the probability to 
add bosonic sources, and on the other by the probability to shift one of the sources and 
hence move within . Because the ‘put/remove’ update step does not alter the fermionic 
bond configuration, we have Zq = Q^{0). On the other hand it adds or removes a fermionic 
monomer term M{(p) at site x. The relative weight of the configurations with or without this 
term is given by Qi/Qo and the acceptance ratios on a lattice with Lt sites become 


Prm(x —>■ 0 ) 


Tput(0 ^ x) 


2 Qo(N(x)) 
Lt Qi(iV(x))’ 

Lt Qi(N(x)) 

2 Qo(N(x)) • 


(14) 

(15) 


The factor Lt compensates for the proposition probability to choose lattice site x out of Lt 
possibilities when putting the sources, while the factor 2 compensates for the asymmetric 
shift proposal probability when moving ip from x to x — 1, since the shift of ip from x to x + 1 
is not allowed. 
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Next we consider the shift update step for the case when the source V' at site x + 1 is 
shifted backwards to site x which is already occupied by the sink ip. The step is graphically 
illustrated in the upper half of figure 4. While the resulting fermion bond configuration is a 

ijj Ip 





Figure 4: Fermionic bond configuration update algorithm. Graphical representation of the hybrid 
‘shift/remove’ update step x + l^x—>-0m backward direction, balanced with the ‘put/shift’ update 
step, 0 ^ X —>■ X + 1. The bosonic background bond configuration is not drawn. 

valid one (it belongs to Zi), the whole fermion conhguration including the source and the sink 
represents an unphysical situation, and in fact does not contribute to any physical observable, 
as discussed before. Therefore, such a backward shift from to Zi, essentially closing the 
open fermion string, automatically induces the removal of the fermionic source and sink pair 
ipip from site x as illustrated in the lower half of figure 4. Such a step is called a hybrid 
‘shift/remove’ update step. Of course, the step is balanced with a hybrid ‘put/shift’ update 
step when the additional fermionic sink and source variables are put on a closed fermion loop 
at the site x. As usual, the acceptance ratios for the hybrid update steps can be read off from 
the weights of the configurations involved and yield 


2 

-Psh/rm(^ + 1^X^0) = —, 

(16) 

-^put/sh (^ ^ ^ ^ 37 -|- 1) 2 * 

(17) 


The factor Lt compensates for the proposition probability to choose the same lattice site x 
when putting the sources ipip back on the lattice, whereas the factor 2 compensates for the 
proposition probability to shift in forward or backward direction when the fermion string is 
still open. Note that there are no ratios of Q-weights involved, since no monomer term is 
added or removed by the hybrid shift/remove update step. 

To complete our discussion of the fermionic bond update, we note that the algorithm 
provides improved estimators for the fermionic two-point function (x) and the partition 
functions Zp. Because the algorithm samples directly the conhguration space Q-^, every open 
fermion string conhguration contributes unity to the stochastic Monte Carlo estimator for 
g^{x). To be precise we have 

g^{x\C G Q^) = 6 ^,xi-x 2 , (18) 

where xi and X 2 are the end and starting point of the open fermion string, i.e., the positions 
of the sink ip and the source ip, respectively. Similarly, every bond conhguration in Zp is 
generated with its proper weight and hence contributes unity to the stochastic estimator for 
Zp, i.e., the Monte Carlo estimator for Zp is simply 

Zp{CeZp) = l. (19) 
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Finally we note that the factors of Lt appearing in the acceptance ratios above may become 
inconvenient in practice, especially towards the continuum limit when —)• oo. The factors 
only occur when contact between Zp and is made, i.e., they are responsible for getting the 
relative normalisation between Zp and right. However, since we make use of translational 
invariance in eq.(18) the factors of Lt are in fact cancelled and can hence be omitted. 


2.2 Updating the bosonic bond configuration 


In this section we now discuss the update steps which relate the bond configuration spaces 
Zp and Q^p for a fixed fermionic bond configuration with fermion number F = 0,1. 

We point out that for an arbitrary superpotential there are in general no restrictions 
on the bosonic bond configurations. This is for example the case for the superpotential Pb 
which we consider in our series of papers, cf. eq.(53). In contrast, the superpotential Pu in 
eq.(54) yields the local constraint N{x) = 0 mod 2 on the site occupation number, due to 
the parity symmetry (j) —>■ —cj). In the following discussion, we always present the generic case 
first, and then specify the modifications or simplifications due to the constraint. In analogy 
to the fermionic bond update, the ‘put/remove’ and the ‘shift’ updates are the main steps 
for updating the bosonic bond configurations. The ‘put/remove’ step introduces or removes 
one or two sources (/>, while the ‘shift’ step shifts the sources by one lattice spacing. If there 
are no restrictions on the bond configuration, we are free to decide for each Monte Carlo 
step whether to proceed by a ‘remove’ update or a ‘shift’ update. With probability prm, we 
propose to remove the sources from the lattice, while the proposition to continue the worm 
update with a ‘shift’ step is chosen with probability 1 — prm- 

The step from Qp to Zp (and vica versa) is induced by removing (or introducing) a bosonic 
source (f) at sites xi and X 2 , with xi = X 2 not excluded. The step does not alter the bond 
configuration, but only the site occupation numbers at sites xi and X 2 - Thus, only the ratios 
of the site weights Qp are involved in the acceptance probability for e.g. the ‘remove’ step, 


Pvm{xi^ X2 y 0) 


' 1 Qp(N(xi)-2) 

PrmLt Qp{N{xi)) 

1 QpjNjxi) - 1) Qp{N{x 2 ) - 1) 
, Qf{N{xi)) Qp{N{x2)) 


( 20 ) 


The prefactor l/(prmTf) is motivated as follows. The factor 1/Tf balances the probability 
for the proposition of putting the bosonic sources at the sites xi and X 2 when re-entering the 
configuration space Q^p, while the factor 1/prm balances the proposition probability for the 
choice of proceeding by the shift update instead of the remove update, as discussed above. 
The acceptance ratios for re-entering the configuration space G^p from Zp are given by 


T’put(0 Xi,X2) 


PrraL 


2 

t 


Qp{N{xi) + 2) 
Qp{N{xi)) 




^ Qp{N{xi) + l)Qp{N{x 2 ) + l) 
* Qp{N{xi)) Qf{N{x 2 )) 


if Xi = X2, 


if Xi / X2. 


( 21 ) 


Two remarks are in order. Firstly, if there are no constraints on the bond configuration, 
one can in principle introduce just a single source cj) which subsequently is shifted around. 
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In effect, the algorithm then samples the one-point function which in this situation is indeed 
nonvanishing. Secondly, we note that if the constraint = 0 mod 2 is in place, the two 
sources can only be placed or removed when xi = X 2 - As a consequence, only the first of the 
two acceptance ratios in eq.(20) and eq.(21) are relevant, while the second ones are zero by 
definition. 

Next, we discuss the bosonic ‘shift’ update. With this step we now change the bosonic 
bond configuration. Shifting the source from site x to a next neighbouring site y is always 
associated with an increase or a decrease of the bosonic bond occupation number between 
the sites x and y by one. Whether or not the occupation number is increased or decreased 
is decided with probability 1/2. Similarly, the source can move forward or backward, and we 
propose both directions with equal probability 1/2. In addition, when there are several types 
of bosonic bonds bi with f G {j —)• k\j,k G N}, we need to decide in each step which bond 
is updated. We do so by choosing the proposition probabilities pj^k with YljkPj^k = 1- 
However, because the proposals are completely symmetric, these probabilities do not affect 
the acceptance ratios. In the following, we will use the shorthand notation 

^ f ^Uk(^) if 2/= 3^ + 1, 

1 rfj^kiy) if 2/ = a: - 1, 

for the occupation number of the bosonic bonds bj^k between the sites x and y. The shifts 
X —)• y and —)> + 1 are balanced with shifts y ^ x and —)> — 1, which 

gives the acceptance ratios 


( 22 ) 


^ y, + 1) 

Wj^k QFi^ix)+j — l) Q_p(A^(y) + A;-|-1) -f _ i -i 

+1 <3f(iV{i)) QHmyf) ' 

1(1, QF(N[x)+k-l) QF[kk(y) +j-FI) _ 

+ 1 Qf(N{x)) ' QF{N(y)) ' 


Tsh(3^ 


y,n 


j^k 

xy 


—)■ n 


j^k 

xy 


1 ) 


/ j^k 
'klxy 

Qf{N[x) - j - 1) 

QF{N{y)-k + l) 

'^j — 

Qf{N{x)) 

QF{N{y)) 

j^k 

'klxy 

QF{N{x)-k-l) 

Qp{N{y)-j + l) 

'' '^3 — 

Qf{N{x)) 

QF{N{y)) 


if y = X + 1, 


if y = X — 1. 


(24) 


Of course, these generic ratios simplify considerably for the specific bonds 6i,iG{l—>'l,l—>' 
2,1 —7- 3} relevant for the superpotentials considered in our series of papers. For example, 
the acceptance ratios for updating the bond 6i_s.i read 


Psh{x y, ^ nl~^^ + 1) 
Psh{x ^ y, - 1) 


wi^i QpjNjy) + 2) 

ni^i + r QF{N{y)) 

<y^ _ Qf{N{x) - 2) 
'wi^i Qf{N{x)) 


(25) 

(26) 


Because the bond is symmetric, there is no need to distinguish whether y = x + lory = x — 1. 
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To complete the discussion of the bosonic bond update, we point out that the algorithm 
again provides improved estimators for the bosonic two-point function gp{x) and the partition 
functions Zp. As in the fermionic case, the algorithm samples directly the configuration space 
Qp with the correct weighting when the sources are present. Therefore, every configuration 
contributes unity to the stochastic Monte Carlo estimator for gp{x), and we have 

gpi^l^F) — ^xi-X2,x , (27) 

where xi and X 2 are the positions of the two sources. Whenever the bosonic update decides 
to remove the sources, we have a configuration m. Zp and hence a contribution of unity to 
the stochastic estimator for Zp, that is, we have 

Zp{CeZp) = l. (28) 

In complete analogy to the fermionic update we note that the factors of Lt appearing in 
the acceptance ratios of the ‘put/remove’ step can be compensated by adjusting the overall 
normalisation of the two-point function, e.g. by making use of translational invariance. 


3 Calculation of the site weight ratios 

In order to calculate the weight of a bond configuration, it is necessary to know the site 
weights 

/ OO 

# (29) 

-OO 

where V((/>) and M (^) depend on the superpotential and the discretisation employed, and 
F = 0,1 is the fermion number, for arbitrary values of the site occupation number n. The 
values of n required in practice are usually limited to O(IO^). However, it turns out that 
even for moderate values of n of order 0(100) the site weights Qp{n) can quickly grow larger 
than 10^^® or more. As a consequence, the calculation of the site weights quickly becomes 
numerically unstable for growing n. In fact, even for simple potentials when the weights 
can be calculated analytically in terms of confluent hypergeometric functions, the numerical 
evaluation of these functions is difficult for large n, and even specialised libraries such as the 
ones available in Wolfram’s Mathematica [22] appear not to be accurate enough. 

Fortunately, for the Monte Carlo simulations we only need ratios of the site weights, such 
as Qp{n + ^)/QFin),Qpin + 2)/Qp{n) and Qi{n) /QQ{n), and these ratios usually do not 
become larger than 0(10) even for large n. In addition, also the transfer matrix elements can 
be rewritten in terms of these ratios as discussed in the appendix of our second paper of the 
series [5]. Therefore, we now present a numerically stable computational strategy to calculate 
the site weight ratios reliably for arbitrary values of the site occupation numbers. 

We start by dehning an arbitrary polynomial superpotential 

p 

Picj)) = (30) 

i=0 

and the corresponding bosonic self-interaction potential V {4>) as well as the monomer weight 

2(p-l) p-2 

¥{(!>)= Y. hcj)\ M(<^) = J^m#. (31) 

^=0 i=0 
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Explicitly, the weights in each sector are then given by 


/ OO 

del) 

-OO 


( 32 ) 


and 


p-2 

Qo{n) = '^miQi{n + i). 


i=0 

For convenience we also define the ratios of the site weights Qpin), 


Rpiji) 

Rpin) 

Rm{n) 


QF{n 1 ) 
Qrin) 
QF{n 2) 
Qpin) 
Qojn) 
Qi{n) 


(33) 


(34) 

(35) 

(36) 


which are used for the acceptance ratios in the Monte Carlo simulations. In principle, only 
the ratios R'i{n) need to be calculated since all other ratios can be derived from those. For 
example, Ri{n) can be expressed in terms of R'i{n) as 


i?i(n) = R'i{n + l)R!i{n ), 


(37) 


but since in some cases (5i(nodd) = 0 the introduction of Ri{2n) is nevertheless necessary. 
Rm{n) can be expressed via the ratios Ri{n) and 7?i(n) and appropriate products thereof, 

Rm{n) = mo+i?i(n) (mi + Ri{n + 2) (m 3 + .. .))+i?i(n) (m 2 + iii(n + 2) (m 4 + ...)), (38) 

and the ratios 7?o(n) and Ro(n) via Rm(n), Ri(n) and R[(n) by 

D f \ 7?m(^ + 2 ) eon\ 

Ro{n) = Ri{n), (39) 

Rm{n) 

r>/e \ Rmin1) / . , /'/in'i 

RoH = , . Ri{n). (40) 


First, we now discuss how to gain numerical stability for the special case of an even 
superpotential P{(!))■ In a second step we will then adapt the idea to treat the somewhat 
more subtle case of an arbitrary superpotential. 


3.1 Even superpotential 

Unbroken supersymmetric quantum mechanics requires a superpotential P{<f)) with deg{P{(p)) 
= 0 mod 2. In particular, in our series of papers we investigate the superpotential 

P{(l)) = P2(p^ (41) 

which is symmetric w.r.t. the parity transformation (j) —)• —(f). As a consequence of the 
symmetry, (^^(f^odd) = 0 for both F = 0,1 and the ratios Rp{n) need not be considered ~ 
instead, it is sufficient to determine Ri{2n) with n G Nq only. 
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For the the potential V ((/>) we then have the form 

, ( 42 ) 

consistent with both the standard discretisation and the Q-exact one. To keep the integrals 
numerically under control, for fixed n we apply a variable transformation (/>—>■(/>/(/) to obtain 
rescaled weights Qi(2n) as 

Qi(2n) = c/)^’^+^Qi(2n). (43) 

Since we have Qi(2n) > 0, we can choose the rescaling factor to be (/> = Qi(2n)^/^^”'’'^) and 
the rescaled weight becomes Qi(2n) = 1. Calculating the ratio of rescaled weights as 

i?,(2n) = = 4(2n + 2), (44) 

Qi(2n) 

where both integrals Qi(2n+2) and Qi(2n) are rescaled with the same factor cj) = Qi(2nY/^‘^^~^^\ 
we find that 

i?i(2n) = ^ i?i(2n). (45) 

In addition, the rescaled weight Qi{2n + 2) is now of 0(1) and can be evaluated reliably 
via numerical integration. So if we start by integrating directly the numerically stable site 
weights Qi(0) and Qi(2), we can recursively generate ratios Ri{2n) with higher and higher 
n. N(^e that after each calculation of a ratio Rii2n), one needs to update the rescaling factor 
(p ^ . This can be achieved most easily via 

~ ~ 2n+l 1 

(j) = (^ 2n+3 _R^(2n) 2"+3 . (46) 

Our procedure guarantees that all involved quantities are of 0(1). Once all ratios Ri{2n) are 
known, one can calculate the ratios Rm{‘^n), noting that for the specific superpotential we 
consider, eq.(38) simplifies to 


Rm{‘2n) = mo + m2Ri{2n). (47) 

The calculation of the ratios i?o(2n) as given in eq.(39) is then straightforward. 

3.2 Arbitrary Superpotential 

In the context of broken supersymmetric quantum mechanics, one encounters superpotentials 
with deg{P{4>)) = 1 mod 2. Therefore, we now adapt the procedure from above to super¬ 
potentials of this form. For simplicity, we restrict ourselves to the odd superpotential we 
consider as the example in our series of papers, 

3 

P{(P) = Y.PicP\ (48) 

i=\ 

If at least one of the coefficients pi and p 2 is nonzero, which is always the case for the 
superpotentials we use, V{(p) reads 

V {(p) = ki(p + k2(p^ + kz(p^ + kicp^, (49) 
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and at least one of the coefficients ki and is nonzero either. This has a two important 
consequences. Firstly, the moments defined in eq.(32) are nonzero for n odd, from which it 
follows that the ratios Rp{n) defined in eq.(36) have to be calculated as well. Secondly, the 
weights Qi(n) are no longer necessarily positive. It turns out, however, that for all practical 
purposes it does not affect the simulations. We will discuss this further in Section 4. 

For the evaluation of the integrals, we apply the same variable transformation cj) ^ <j)/^ 
as before, such that we have rescaled weights Qi{n) given by 

Qi(n) = 0”+^ Qi(n). (50) 

We now choose cj) = ■ sgn((5i(n)). Then, the integral becomes Qi(n) = 1 again 

as before. Furthermore, defining the rescaled ratios R!i{n) to be 

^((n) = = Qi(n +1), (51) 

Qi{n) 

where both integrals Qi(n+1) and Qi{n) are rescaled with the same factor (j) = |(5i(n)|^/("'+^) • 
sgn(Qi(n)), we find R'i{n) = (t)R\{n). We proceed analogously to the case of the even 
superpotential by recursive iteration, with the only exception that we generate the ratios 
R'i{n) instead of the ratios i?i(n). The update for the rescaling factor (/>—)•(/)' is done via 

4> = |(/)| "+2 |i2(^(n)|’"+2 • sgn(i?^(n)). (52) 

Once all the ratios R'i{n) are known, one can calculate the ratios Ri{n) via eq.(37), the ratios 
Rm{n) via eq.(38), and the ratios Ro{n) and R^in) via eq.(39) and (40), respectively. 

4 Results of the Monte Carlo simulations 

The results in this section are merely thought of as a proof of the feasibility of the algorithm 
and as a test of its efficiency. Comparing the Monte Carlo results with the exact solution of 
the system at finite lattice spacing provided in our second paper [5] of course also serves as 
a validation for the algorithm. We refer to that paper for a thorough discussion and physical 
interpretation of the results. 

For the following Monte Carlo simulations, we consider the same superpotentials and 
discretisations as in the previous two papers. In particular, we simulate the system using 
the action with counterterm for both unbroken and broken supersymmetry as well as the 
Q-exact action for unbroken supersymmetry. Details for the various actions can be found in 
the first paper of our series. Here we only give the details of the superpotentials for broken 
and unbroken supersymmetry, respectively, 

n(<?^) = -^</>+^A,/.3, (53) 

, (54) 

and we recall that the continuum limit is taken by fixing the dimensionful parameters /U, g, A 
and L while taking the lattice spacing a —)• 0. In practice, the dimensionless ratios fu = 9 / 
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Figure 5: Unbroken supersymmetric quantum mechanics, standard discretisation. Bosonic correlation func¬ 
tion for antiperiodic (black) and periodic b.c. (red) (lying on top of each other) and fermionic correlation 
function for antiperiodic (green) and periodic b.c. (blue) for fj,L = 10 at coupling /„ = 1. The dashed lines are 
the exact results from [5]. 


and fb = fix the couplings and /tL the extent of the system in units of while 

and a/L are subsequently sent to zero. In analogy to the number of sweeps for a standard 
Monte Carlo simulation, we count the number of times the algorithm is in either one of the 
two configuration spaces 2f,F = 0,1. The statistics for a simulation is therefore given by 

■^0 + -^1 = 2 a . 

First, we consider the standard discretisation with the superpotential Pu such that super- 
symmetry is unbroken. As a first observable, we show the results for the bosonic and fermionic 
correlation functions for fiL = 10, L/a = 60 and fu = ^ for Za = 10^ in figure 5. This is 
essentially the same plot as figure 10(b) in our second paper [5], but now with the additional 
data from the Monte Carlo simulation and plotted on a logarithmic scale. Note that we use 
the notation x = t in accordance with [5] . The simulation indeed reproduces the exact result 
within very small statistical errors which demonstrate the efficiency of the algorithm. The 
exponential error reduction is due to the use of the improved estimators for the two-point 
function which are available in the context of the worm algorithms. The improvement is par¬ 
ticularly impressive for the fermionic correlator where the error reduction allows to follow the 
correlator over more than seven orders of magnitude without loss of statistical significance. 
In fact the relative error for the lowest value of the fermionic correlator is still only 4%. 

As a second example, we show the mass gaps for different at a coupling fu = ^ with 
statistics of Za = 10® in figure 6. The /uL considered are in the region where thermal effects 
are negligible and essentially only Zq contributes to the total partition function, such that 
Za — ZQ. We extract the masses from the asymptotic behaviour of the correlation function at 
large t, i.e., we extract the lowest energy gap. Because of the extremely good signal-to-noise 
ratio the asymptotic behaviour can be truly reached and, in doing so, systematic errors from 
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Figure 6: Unbroken supersymmetric quantum mechanics, standard discretisation. Continuum limit of the 
lowest bosonic (circles) and fermionic (squares) mass gap for fiL — 10 (black), jiL — 19 (red), pL = 31 (blue) 
and /„ = 1. The inset shows a detailed comparison with the exact results (dashed lines). 


contributions of excited states are essentially excluded. Of course, we know from our exact 
results that the overlap of the simple operators we use to construct the two-point function 
is close to maximal. This is clearly visible in figure 5 where we observe an almost purely 
exponential decay for all t/L. Because the energy gaps are independent of /rL, they are 
expected to fall on top of each other for all values of fiL at fixed lattice spacing o/i. This is 
indeed the case within our numerical accuracy, and the extracted masses, when expressed in 
units of /U, indeed extrapolate to the correct zero-temperature continuum limit. The inset of 
figure 6 shows a detailed comparison of the simulation results with our exact solution from 
[5] represented by the dashed line and we observe a beautiful agreement even very close to 
the continuum. 

Next, we consider the action with counterterm and the superpotential Pf, for which the 
supersymmetry is broken. In this case we encounter an issue concerning the potential non¬ 
positivity of the weights which we already mentioned in Section 3.2. This potentially danger¬ 
ous sign problem is not of fermionic origin, but is instead related to the bond formulation of 
the bosonic degrees of freedom. As a matter of fact it occurs already in the purely bosonic 
system, independent of the dimensionality of the system. However, negative weights only 
occur in a region of parameter space which becomes irrelevant towards the continuum limit. 
In that sense, the sign problem is a lattice artefact and can be avoided straightforwardly. 
Nevertheless, in order to eliminate any systematic error we deal with this bosonic sign prob¬ 
lem by incorporating the sign of the configuration into the observables, even though it has no 
practical consequences. 

As a first observable in the broken case, we show the bosonic and fermionic two-point 
functions, {(j)t4>o) and for periodic and antiperiodic b.c. for fiL = 10 at fixed coupling 
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Figure 7: Broken supersymmetric quantum mechanics, standard discretisation. The bosonic two-point func¬ 
tion for antiperiodic (black) and periodic b.c. (red) and the fermionic one for antiperiodic (green) and periodic 
b.c. (blue) for L/a = 60, /iL = 10 at coupling fi, = 1. The dashed lines are the exact results from [5]. 

/fe = 1 in figure 7 for a statistics of Za = 10®. The exact results from [5] are shown as dashed 
lines. The simulation yields results which agree with the exact results within the very small 
statistical errors on the level of l%o. Note that the correlators for periodic and antiperiodic 
b.c. are constructed a posteriori from the simulation results in the bosonic and fermionic 
sectors Zq and Zi, respectively, and it is crucial to sample the relative weight between the 
two sectors correctly in order to get the final values right. The relative sampling is solely in 
the responsibility of the fermion simulation algorithm. Our results in figure 7 show that the 
open fermion string algorithm indeed transits sufficiently well between the two sectors. 

This statement can be made more quantitative by looking at the ratio ZpjZa which rep¬ 
resents the Witten index in our field theoretic setup. From our exact results in [5] we expect 
a nonzero Witten index at finite lattice spacing which however extrapolates to zero in the 
continuum limit. So the behaviour of the algorithm towards the continuum limit is particu¬ 
larly interesting, because for vanishing lattice spacing the would-be Goldstino at finite lattice 
spacing turns into a true, massless Goldstino. In such a situation one usually encounters 
critical slowing down of the simulation algorithms, such that the errors on the results grow 
large and the results become unreliable. The massless Goldstino is directly related to the tun¬ 
nelling between the bosonic and the fermionic sector, and the reproduction of a Witten index 
VF = 0 in the continuum with small errors is hence a true demonstration of the efficiency of 
the open fermion string algorithm to transit between the bosonic and fermionic sector. In 
addition, we know from [5] that the lattice artefacts are exponentially enhanced towards zero 
temperature and it is interesting to see how the simulation algorithm handles this situation 
at coarse lattice spacing. 

In figure 8 we show the ratio Z^jZa as a function of the lattice spacing afj, for different 
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Figure 8: Broken supersymmetric quantum mechanics, standard discretisation. Continuum limit of the 
partition function ratio Zp/Za, i.e., the Witten index, for /rL = 5 (blue), jiL — 10 (red), ^iL = 20 (green), 
jiL = 30 (black) at fixed coupling fb — 1. The dashed lines are the exact results from [5]. 


values of fiL at fixed coupling /{, = !. For this quantity, too, the simulation yields results 
which agree with the exact results within the small statistical errors. Moreover, the efficiency 
of the algorithm does not appear to deteriorate towards the continuum limit or for small 
values of fj,L where the Witten index is very close to zero. This can for example be seen from 
the fact that the errors obtained with fixed statistics essentially remain constant towards the 
continuum limit and are also independent of the system size. This nicely demonstrates the 
efficiency of the algorithm also for a system with broken supersymmetry. 

The last system we investigate with the worm algorithm is unbroken supersymmetry 
formulated with the Q-exact action^. We first consider the ratio of partition functions ZpjZa 
which in the limit of [iL —?■ oo yields the Witten index. From a simulational point of view, 
the ratio essentially calculates the fraction of configurations in sector Zq versus the ones in 
Z\. For unbroken supersymmetry the system is almost exclusively in the bosonic sector, and 
hence the ratio is very close to one except when the size of the system becomes very small, i.e., 
in the high temperature limit. Moreover, from our exact results in [5] we know that the lattice 
artefacts in this quantity are very small and the continuum limit is not very interesting. For 
these reasons, we consider in figure 9 the dependence of the ratio Z^jZa on [xL for different 
values of the lattice spacing ajL with a statistics of Za = 10®. 

Also for this quantity, we find that the results agree with the exact result within the very 
small statistical errors. Again, the open fermion string algorithm proves to be very efficient 
even close to [xL ~ 0 where the tunnelling from the bosonic to the fermionic sector and vice 

^For Monte Carlo simulations using the Q-exact action for broken supersymmetry, we encounter the very 
same problems we ran into in the transfer matrix approach. The bond occupation number grows extremely 
large even on small lattices and for coarse lattice spacings such that the generation of reliable results turns out 
to be impossible. 
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Figure 9: Unbroken supersymmetric quantum mechanics, Q-exact discretisation. Z^ijZa as a function of 
for Lja — 16 (red), L/a = 32 (blue), L/a = 64 (black) at fixed coupling /„ = 1. 


versa becomes important and dominates the behavionr of the system. Thus, even in this 
somewhat extreme situation of very high temperature, the algorithm does not show any signs 
of critical slowing down despite the fact that there is a quasi-zero mode in the system. 

Note that the algorithm is capable of handling negative bare masses independent of the 
discretisation used and fig. 9 is simply also an illustration of this fact. 

The last quantity we calculate are the lowest bosonic and fermionic mass gaps for different 
at fixed coupling of /„ = 1 from a statistics of Za = 10®. The mass gaps are extracted 
from the two-point correlation functions exactly in the same way as before for the standard 
action, and in figure 10 we show the results of this analysis. As expected, the masses for the 
boson and the fermion are indeed indistinguishable within statistical errors. The degeneracy 
of the masses at hnite lattice spacing due to the Q-exactness of the action emerges also for the 
results from Monte Carlo simulations. Note that the chosen values for lie well within the 
region where thermal effects are negligible and the masses extrapolate nicely to the correct 
zero-temperature continuum limit. The inset in figure 10 shows a detailed comparison with 
our exact results from [5] and we again observe beautiful agreement. 


5 Conclusions 

In this paper we present an algorithm for simulating M = 2 supersymmetric quantum me¬ 
chanics on the lattice. The algorithm is based on the reformulation of the system in terms of 
bosonic and fermionic bonds, and in essence represents an efficient Monte Carlo scheme for 
updating fermionic and bosonic bond configurations. The updating of the fermionic degrees 
of freedom is of specific interest, because this is in general the most challenging part of a 
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Figure 10: Unbroken supersymmetric quantum mechanics, Q-exact discretisation. Continuum limit of the 
lowest bosonic (black squares) and fermionic (red circles) mass gaps for ^.L = 17, and bosonic (blue squares) 
and fermionic (green circles) mass gaps for = 31 at fixed coupling /u = 1. The inset shows a detailed 
comparison with the exact result (dashed line). 


simulation. This is particularly true for systems with broken supersymmetry, where standard 
simulation algorithms suffer from critical slowing down due to the massless Goldstino mode. 
In addition, these systems inevitably also suffer from a sign problem related to the Goldstino 
and the vanishing Witten index. 

In contrast, the fermion simulation algorithm proposed in [2] eliminates critical slowing 
down by directly sampling the fermionic two-point correlation function. It is based on intro¬ 
ducing a fluctuating open fermion string which efficiently updates the bond configurations on 
all length scales up to the correlation length associated with the fermionic correlation function. 
As a consequence, the fermion string induces frequent tunnellings between the bosonic and 
fermionic vacuum when that correlation length becomes large. Since the two vacua contribute 
to the partition function with opposite signs, the frequent tunnelling guarantees sufficiently 
small statistical fluctuations for the average sign, and hence a solution to the fermion sign 
problem. In fact, the more severe the sign problem gets towards the continuum limit, the 
more efficiently the algorithm tunnels between the bosonic and fermionic sectors. This is of 
course due to the growing correlation length associated with the vanishing Goldstino mass. 
The bosonic degrees of freedom can be expressed in terms of bonds as well. Therefore, we 
also give the details of an updating algorithm for the bosonic bond configurations. Since 
we consider Q-exact discretisations in addition to the standard one, the algorithm involves 
updating generic types of bonds. 

The simulation algorithm requires the calculation of the site weights Qf{N). Their nu¬ 
merical evaluation, however, turns out to be numerically unstable for large site occupation 
numbers N. Hence, in Section 3, we devise a computational strategy which allows to reliably 
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evaluate the ratios of weights for arbitrarily large occupation numbers. Since this is a generic 
problem occurring in the bond formulation of field theories with real scalar fields, such a 
computational scheme is useful also in other sitnations. 

Finally, we present a selection of results obtained using the open fermion string algorithm. 
We concentrate on two specihc realisations of supersymmetric quantum mechanics, one with 
broken and one with unbroken supersymmetry. In addition, we consider both the standard 
and the Q-exact discretisation. Since exact results are available at finite lattice spacing 
from our investigation in [5], we can benchmark our stochastic results and directly validate 
them. The calculation of the bosonic and fermionic correlation functions shows that they 
can be determined very accurately over several orders of magnitude. This allows for a very 
precise computation of the boson and fermion masses, the latter in many cases with a smaller 
error than the former. In general, a precision of l%o can be reached with a very modest 
computational effort. In systems with broken supersymmetry it is crucial that the simulation 
algorithm efficiently samples the relative weights between the bosonic and fermionic sectors. 
Our results for the partition function ratio i.e., the Witten index, show that this is 

indeed the case. For fixed statistics, the errors do not grow towards the continuum limit. In 
that limit the index gets very close to zero and the sign problem would therefore be most 
severe. Similarly, the error is essentially independent of the system size, which shows that 
the sign problem is truly solved. 
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